rm(list = ls())
gc()

library(lubridate)
library(foreign)
library(ggplot2)
library(ggthemes)
library(estimatr)
library(survey)
library(tidyverse)
library(mice)
df = readRDS("data/baseline_clean_mi.rds")

input = svydesign(ids = ~1, probs = NULL, weights = df$weights, data = df)

df$rtrn_head_2years[df$rtrn_head_2years < 1] %>% table
df$rtrn_head_12[df$rtrn_head_12 < 1] %>% table
df$rtrn_head_ever %>% table(useNA = "always")

m = 10
n = 3003

# Figure 1: descriptive--------------------------------------------------------------------
df_plot = df %>% select(response_num, .imp, rtrn_head_12, expect2yr_syria, rtrn_head_ever, weights)

df_plot = df_plot %>% 
  pivot_longer(-c(response_num, weights, .imp))

df_plot = data.frame(df_plot)
vars = unique(df_plot$name)
store = matrix(NA, ncol = 3, nrow=length(vars))
colnames(store) = c("var", "mean", "sd")
rownames(store) = vars

means = NULL
variances = NULL
for(i in 1:length(vars)){
  for(j in 1:m){
    df_temp = df_plot[df_plot$name == vars[i] & df_plot$.imp == j,]
    des = svydesign(id=~1, weights=~weights, data = df_temp)
    means[j] = svymean(~value, design = des)
    variances[j] = svyvar(~value, design = des)
  }
  store[i, 1] = vars[i]
  store[i, 2] = pool.scalar(means, variances, n = n, k = 1)$qbar  # Barnard-Rubin 1999. Same as mean(means)
  store[i, 3] = pool.scalar(means, variances, n = n, k = 1)$t %>% sqrt  # Barnard-Rubin 1999. Same as within variance [mean(vars)] + between variance [sum((means - mean(means))^2)/(m - 1)] + between variance/m [sum((means - mean(means))^2)/(m - 1)/m]
}

df_plot2 = data.frame(variable= rownames(store),store[,2:3]) %>% 
  mutate(mean = mean %>% as.character() %>% as.numeric(),
         sd = sd %>% as.character() %>% as.numeric())

df_plot2$low = df_plot2$mean - 1.96*df_plot2$sd/sqrt(n)
df_plot2$up = df_plot2$mean + 1.96*df_plot2$sd/sqrt(n)

df_plot2 = df_plot2 %>% 
  mutate(variable2 = c("Return in 1 year", "Return in 2 years", "Return ever"),
         variable2 = factor(variable2, levels = c("Return in 1 year", "Return in 2 years", "Return ever")))

df_plot2 %>% 
  ggplot(aes(x = variable2, y = mean, ymin = low, ymax = up, fill = "a")) + 
  geom_col() + 
  geom_errorbar(width = .3) + 
  theme_bw() + 
  scale_fill_hc() + 
  theme(legend.position = "none") + 
  labs(x = "", y = "") + 
  scale_y_continuous(labels = scales::percent) + 
  scale_x_discrete() + 
  annotate(geom = "text", x = df_plot2$variable2, y = df_plot2$mean + .04, label = scales::percent(df_plot2$mean))
ggsave("figures/return_intentions.png", height = 4, width = 4)
